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We make simulations with 2 flavor Wilson fermions to investigate the nature of the end points of 
Roberge- Weiss (RW) first order phase transition lines. The simulations are carried out at 9 values 
of the hopping parameter k ranging from 0.155 to 0.198 on different lattice spatial volume. The 
Binder cumulants, susceptibilities and reweighted distributions of the imaginary part of Polyakov 
loop are employed to determine the nature of the end points of RW transition lines. The simulations 
show that the RW end points are of first order at the values of k in our simulations. 

PACS numbers: 12.38.Gc, ll.10.Wx, 11.15.Ha, 12.38.Mh 



I. INTRODUCTION 



A full understanding of QCD phase diagram is of great 
importance theoretically and phenomenologically. The 
QCD phase diagram addresses which forms of nuclear 
matter exist at different finite temperature and baryon 
density, and whether there are bona fide phase transi- 
tion separating them, thus it is essential for relativistic 
heavy ion collision experiments and astrophysics. QCD 
is a strongly interacting theory on the scales of a baryon 
mass and below, so non-perturbative calculations from 
the first principle are preferrable. Despite that substan- 
tial progress has been made with Monte Carlo simula- 
tions of lattice QCD at zero baryon density, the stud- 
ies at nonzero baryon density are haunted by the "sign" 
problem, for example, see Ref. [l|. To date many indi- 
rect methods have been proposed to circumvent the sign 
problem, overviews with references to these methods can 
be found in Ref. [I], One of these methods consists 
of simulating QCD with the imaginary chemical poten- 
tial for which the fermion determinant is positive |3l— KLlj) - 
Full information can be obtained by using the imaginary 
chemical potential which allows for analytic continuation 
via truncated polynomials. 

The phase structure of QCD with imaginary chemical 
potential not only deserves detailed investigations in its 
own right theoretically, but also has significant relevance 
to phy sics at zero or small real chemical potential 0-0, 
ll2Hl6| . QCD with imaginary chemical potential has a 
rich phase diagram as a function of imaginary chemical 
potential and quark masses. 

In this paper, we present a study of phase structure of 
QCD at fixed imaginary chemical potential 9 = /T = 
7r for Nf = 2 QCD with Wilson quarks. The partition 



function including the imaginary chemical potential is 
Z(T,p J )=Tr(e-M H -* li '* r >\ (l) 



Roberge and Weiss made the essential work with the 
imaginary chemical potential (l7j . they found that the 
partition function of QCD with imaginary chemical po- 
tential has two important symmetries, reflection symme- 
try in [i — /ifl -H/i/ and periodicity in imaginary chemical 
potential, 



Z(T,fi) = Z(T ) - li ), 
Z(p/T) = Z{n/T + i2nn/3). 



(2) 
(3) 
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The periodicity is smoothly realized in the low temper- 
ature, strong-coupling regime, whereas in the high tem- 
perature, weak-coupling regime, it is realized in a non- 
analytic way. At high temperature, the system undergoes 
a first order transition (RW transition) at critical values 
of the imaginary chemical potential fii /T = (n + ^)2w/3 
[l7l - [l9j between adjacent Z(3) sectors and these Z(3) sec- 
tors are characterized by the Polyakov loop. Thus the 
picture for the T — 9 phase diagram is that repeats with a 
periodicity the first order transition line in the high tem- 
perature regime which necessarily ends at an end point 
at some temperature Trw when the temperature is de- 
creased sufficiently. 

At these end points, there are evidence that the ana- 
lytic continuation of deconfinement/chiral transition line 
from real chemical potential to imaginary chemical one 
meets the RW transition line. Recent numerical studies 
show that the RW transition line end points are triple 
points for small and heavy quark mass, and second or- 
der end points for intermediate quark masses. So there 
exist two tricritical points which separate the first order 
regime from the second one [3 5]. Moreover, it is pointed 
out p| EH EH that the scaling behaviour at the tricriti- 
cal points may shape the the critical line for real chem- 
ical potential, and subsequently, the line for real chemi- 
cal potential is qualitatively consistent with the scenario 
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suggested in Ref. which show that the first order 

region shrinks with the increasing real chemical potential. 

Most of studies of finite temperature QCD have been 
performed usin g stagg ered fermion action or the im- 
proved versions [2CH26l [ , staggered fermion approach and 
Wilson fermion approach have their own advantages and 
disadvantages, for example, see Ref. [13] . The KS fermion 
formalism preserves the U(l) chiral symmetry, whereas 
it needs a fourth root trick for one flavour which might 
lead to locality problem [28| and phase ambiguities [23 |. 
On the contrary, Wilson fermions completely solve the 
species doubling problem, whereas it suffers from an ex- 
plicit chiral symmetry breaking. The lattice simulation 
with Wilson fermions is more time-consuming than stag- 
gered fermions, it can provide complementary informa- 
tion and crosscheck to simulations with other actions and 
establish a better understanding of QCD phase diagram. 

In this paper, we attempt to investigate the RW transi- 
tion line end points by lattice QCD with two degenerate 
flavors of Wilson fermions. In Sec.|TTl we define the lattice 
action with imaginary chemical potential and the physi- 
cal observables we calculate. Our simulation results are 
presented in Sec. Mil followed by discussions in Sec. [IVJ 



II. LATTICE FORMULATION WITH 
IMAGINARY CHEMICAL POTENTIAL 



We consider the partition function of system with 
Nf = 2 degenerate flavors of Wilson quarks with chemi- 



cal potential on the lattice 



Z = J [dU}[dip}[dip}e- b °~ b f 
= I [dU] [ BetM[U, 



(4) 



where S g is the gauge action, and Sf is the quark action 
with the quark imaginary chemical potential /!/ = 9T. 
For S g , we use the standard one-plaquette action 



(5) 



where f3 = 6/g 2 , and the plaquette variable U p is the 
ordered product of link variables U around an elementary 
plaquette. For Sf, we use the the standard Wilson action 



The fermion matrix is 

3 

M XyV (U,K,fi) = 8 x>y - K } j 

3 = 1 



+(l + 7i)tfj t (*-5X I v+3 



(l- li )e a ^U i (x)8 x 



+{l + li )e^Ul{x-l)5 x 



a+4 



(7) 



We carry out simulations at 9 = ir. As it is pointed out 
that the system is invariant under the charge conjugation 
at 6 = 0,7r, when 9 is fixed [l4j. But the 0-odd quantity 
0(9) is not invariant at 9 = ir under charge conjugation. 
When T < T R w, 0(9) is a smooth function of 9, so it is 
zero at 9 = n. Whereas when T > T R w, the two charge 
violating solutions cross each other at 9 = jr. Thus the 
charge symmetry is spontaneously broken there and the 
0-odd quantity 0(9) can be taken as order parameter . 
In this paper, we take the imaginary part of Polyakov 
loop as the order parameter. 

The Polyakov loop L is defined as the following: 



(L) 



'N t 

l[u 4 ( x ,t) 



(8) 



here and in the following, V is the spatial lattice vol- 
ume. To simplify the notations, we use X to represent 
the imaginary part of Polyakov loop L, X — Im(L). 

The susceptibility of imaginary part of Polyakov loop 
X is defined as 



x = v((x-(x)f) 

which is expected to scale as: d, Q 
X = Llf»cj>{TL 1 J»), 



(9) 



(10) 



where r is the reduced temperature r = (T — Trw)/Trw, 
V = Ll . This means that the curves x/lV" at differ- 
ent lattice volume should collapse with the same curve 
when plotted against tL 1 J v ' . In the following, we em- 
ploy f3 - 13 RW in place of r = (T - T RW )/T RW . The 
critical exponents relevant to our study are collected in 
Table. U[i,|30|]. 





v 7 7/1/ 


3D ising 
tricritical 
first order 


0.6301(4) 1.2372(5) 1.963 
1/2 1 2 
1/3 1 3 



N f 

s f = Y,Y,^ M *>y( u > K ^)My), (6) 

/=! x,y 

where k is the hopping parameter, related to the bare 
quark mass m and lattice spacing a by k = 1/ (2am + 8). 



TABLE I: Critical exponents relevant to our study. 

We also consider the Binder cumulant of the imaginary 
part of Polyakov loop which is defined as the following: 



B 4 = ((X-(X))*)/((X-(X)) 2 Y 



(11) 
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FIG. 1: Scaling behavior of susceptbilities of the imaginary part of the Polyakov loop according to the first order critical indexes 
(the left panel), and to the 3D Ising critical indexes (the right panel) at k = 0.155. 
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FIG. 2: Reweighted distributions of the imaginary part of the Polyakov loop Im(L) at the corresponding end point Prw, an d 
P > Prw and f3 < Prw on each lattice spatial volume at k = 0.155. 



with (X) = 0. In the thermodynamic limit, B^(j3) takes 
on the values 3, 1.5, 1.604, 2 for crossover, first order 
triple point, 3D Ising and tricritical transitions, respec- 
tively. However, on finite spatial volumes, the steps are 
smeared out to continuous functions. In the vicinity 
of the RW transition line end points, B± is a function 
of x = (0 — 0i{w)L l J v and can be expanded as a se- 
ries HEEm. 



B± = B 4 (f3 c , oo ) + aix + a 2 x 2 H , (12) 



III. MC SIMULATION RESULTS 



In this section, we will present our results for sim- 
ulating QCD with two degenerate flavors of Wilson 
fermions at finite temperature T and imaginary chemical 
potential ifii. Both the <fi algorithm with a Metropolis 
accept/reject step and the R algorithm are used [3lj |. 
The simulations are performed on lattice with different 
spatial volume with temporal extent N-t = 4 at k = 
0.155, 0.160, 0.165, 0.168, 0.170, 0.175, 0.180, 0.190, 0.198 
For each k value, we carry out simulations on lattice 
of size L s — 8, 12, 16, and for some n values, lattice of 
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FIG. 3: Scaling behavior of susceptibilities, and Binder cumulants of the imaginary part of the Polyakov loop according to the 
first order critical indexes (left panels), and to the 3D Ising critical indexes (right panels)at k = 0.198. 
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FIG. 4: Reweighted distributions of the imaginary part of the Polyakov loop at k = 0.198 at the corresponding end points 
Prw- 



size L s = 10 or/and L s = 20 are also used. Simulations are carried out with (f> algorithm with a Metropolis 
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accept/reject step with the acceptance rate ranging 
from 42 — 93%, The other simulations are carried out in 
terms of R algorithm with the molecular dynamics time 
step St — 0.01. Ref. [3l[ pointed out that R-algorithm 
has errors of order 0(St 2 ), so the correct results of 
this algorithm consists of extrapolation to zero stepsize. 
However, in practice a short-cut without extrapolation is 
used. Recently, the exact RHMC algorithm is invented 
which also allows many improvements (32| . In our 
simulation, St — 0.01 is sufficiently smaller compared 
with the statistical errors of our simulations. There are 
20 molecular steps for each trajectory. We generate 
20,000 trajectories after 10,000 trajectories as warmup. 
Ten trajectories are carried out between measurements. 
We use the conjugate gradient method to evaluate the 
fermion matrix inversion. 

On each lattice size, we make simulations at typically 
4-6 different ft values. For fixed = iirT, there is tran- 
sition in T between the low temperature phase and the 
high temperature phase. In order to determine the RW 
transition line end point Prw from the peak of suscep- 
tibilities, we use the data obtained through simulations 
at the 4-6 different j3 values, and calculate susceptibili- 
ties at additional /3 values, by employing the Ferrenberg- 
Swendsen reweighting method [33] . 

Let us first present the critical couplings Prw on dif- 
ferent spatial volume at different k in Table. HU 

TABLE II: Results of critical couplings Prw on different spa- 
tial volume at different k, we also make simulations on lattice 
8 3 x 4 at k = 0.185, 0.195, the critical couplings Prw are 
4.8810(20), 4.6610(20), respectively. 



K 


8 


10 


12 


16 


20 


0.155 


5.4319(40) 




5.3887(40) 


5.427(10) 


5.4289(50) 


0.160 


5.361(50) 


5.365(30) 


5.347(10) 


5.3499(60) 




0.165 


5.2566(90) 5.262(13) 5.2493(20) 


5.2412(10) 


5.2581(10) 


0.168 


5.206(15) 




5.2103(22) 


5.2167(6) 


5.2181(10) 


0.170 


5.1645(50) 




5.1722(10) 


5.1770(5) 


5.1785(2) 


0.175 


5.0781(30) 




5.0838(50) 


5.0882(40) 


5.1095(30) 


0.180 


4.9802(20) 




5.0388(60) 


5.0391(40) 




0.190 


4.7800(20) 




4.7658(10) 


4.7883(3) 




0.198 


4.5910(20) 




4.5955(10) 


4.5980(2) 





The presence of a first order phase transition at the 
end point of Roberge- Weiss transition line at K = 0.155 
can be found from the scaling behavior of the susceptibil- 
ities of the imaginary part of Polyakov loop x presented 
in Fig. [T] From Fig. Q] we can find that the rescaling 
quantities x/Ll plotted against (j3 — Prw)L\' v does 
not fall on the same curve completely, whereas peaks of 
the rescaling quantities x/lV" obviously exhibit scal- 
ing behavior which conforms to the first order transition. 
From Eq. (|10p . we can find that the index 7/1/ regulates 
the height of peaks while the index v regulates the width 



of peaks. As a comparison, we also present the behavior 
according to the 3D Ising transition index in the right 
panel of Fig. [T] from which we can find that large devia- 
tion from the 3D Ising scaling behavior manifest clearly. 
At k = 0.160, similar observations of susceptibilities as 
those at ft = 0.155 can be found. 

In Fig. [2J we present reweighted distributions of the 
imaginary part of Polyakov loop Im{L) at the corre- 
sponding Prw and two p values on lattice size L s = 
8, 16, 20. On each lattice size, at Prw, reweighted dis- 
tribution of Im(L) exhibits two-state signal, while, at 
P > Prw and P < Prw, reweighted distributions of 
Im(L) do not exhibit two-state signal. At other n val- 
ues, reweighted distributions of the imaginary part of 
Polyakov loop Im(L) at the corresponding Prw, P > 
Prw an d P < Prw on each lattice size have the same 
observations as those at k — 0.155. For clarity, we only 
present the result at k = 0.168 in the following. 

We also make simulations at k = 0.190, 0.198, the re- 
sults of simulations at k — 0.198 are presented in Fig. [31 
and Fig. [U From the two upper panels of Fig. [3J we 
can find that the first order transition indexes are more 
suitable to describe the behavior than the 3D Ising ones. 
This situation can be made clearer when we look at the 
B4 behavior depicted in down panels of Fig. [3J from which 
we can find that the quantities of Binder cumulant plot- 
ted against rescaling P fall on the same curve completely. 
Note that from Eq. (fTTj) . the scaling behavior of Binder 
cumulants is governed by the critical index v which also 
determines the width of peaks of the rescaling quantities 
x/ Ll^ u . The fact that the value of v for first order tran- 
sition accounts for the width of peaks of x/lV" better 
than the second order transition show that the transi- 
tion is first order, and this situation is confirmed by the 
scaling behavior of Binder cumulant B4. We also present 
reweighted distribution of the imaginary part of Polyakov 
loop at Prw at /t = 0.198 in Fig. 2] which exhibits two- 
state signal. At k — 0.190, similar observations as those 
at k = 0.198 can be observed. 

The results of simulations at k = 0.170, 0.175, 0.180 
are shown in Fig. [5] In view of the fact that large finite- 
size corrections are observed in simple spin models even 
when the transition is first order pj. |34| . we can find that 
the first order transition indexes perform much better 
than the second order transition ones. This observation 
can be enhanced from the reweighted distribution of the 
imaginary part of Polyakov loop presented in Fig. O 

Comparing to the above results, it is difficult to de- 
termine the nature of RW transition line end points at 
k = 0.165, 0.168 results of which are presented in Fig. [JJ 
However, when we look at the behavior at large lattice 
size presented in Fig. [7] it is a reasonable conclusion 
that the behavior of RW transition line end points at 
K = 0.165, 0.168 are of first order. This conclusion can 
be enhanced when we look at the reweighted distributions 
of Im{L) at the end point Prw at K = 0.165 presented 
in Fig. [5[ and reweighted distributions of Im(L) at the 
corresponding Prw, P > Prw and P < Prw ° n lattice 
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FIG. 5: Scaling behavior of susceptibilities of the imaginary part of the Polyakov loop according to first order critical indexes 
(left panels) and to the 3D Ising critical indexes (right panels). 



size L s = 12, 16, 20 at k = 0.168 presented in Fig. [9] 



IV. DISCUSSIONS 



We have studied the nature of critical end points of 
Roberge- Weiss transition of two flavor lattice QCD with 
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FIG. 6: Reweighted distributions of the imaginary part of the Polyakov loop at ft = 0.170, 0.175, 0.180 at the corresponding 
end points Prw- 



Wilson fermions. When = iirT, the imaginary part of 
Polayakov loop is the order parameter for studying the 
transition from low temperature phase to high tempera- 
ture one. Within the imaginary chemical potential for- 
mulation, the partition function is periodic in imaginary 
chemical potential. The different Z(3) sectors are char- 
acterized by the phase of Polyakov loop. The Roberge- 
Weiss transition which occurs at f-ii/T = 2n(k + l/2)/3 
is of first order in the high temperature phase, whereas 
it is of crossover in the low temperature phase. 

Our simulations are carried out at 9 values of n on dif- 
ferent 3-4 spatial volumes. Our lattice N t — 4 is coarse. 
In Ref. j35[, the lattice spacing with 2 flavor Wilson 
fermions at /3 = 5.3 is estimated to be 0.12 — 0.13 fm. In 
Ref. [36| , the lattice spacing with 2 flavor Wilson fermions 
is estimated to be 0.246 fm which is found almost inde- 
pendent of /? in the range of /3 = 3.0 — 4.7. In our simu- 
lations, (3 varies roughly from 4.6 to 5.4, thus, the lattice 
spacing a is estimated to be a ~ 0.12 — 0.25 fm. 

In order to estimate the pseudo-scalar meson mass m w , 
the vector meson mass m p and ratios m,/m p , T c /m p at 
our simulation points, we use the data in Table II in 



Ref. [37} ■ By using the standard quark and gauge action, 
Bitar et al. studied hadron thermodynamics with Wil- 
son fermions on lattice 8 3 x 4 and calculated the zero 
temperature hadron mass on lattice 8 3 x 16 with dy- 
namical fermions. We compile their results and present 
in the following: at k = 0.16, (3 — 5.28, m^jra p = 
0.943(3), T c /m p = 0.19425(7), at ft = 0.17, £ = 5.12, 
m^/rrip = 0.899(4), T c /m p = 0.2066(8), at « = 0.18, (3 = 
4.94, m^/iJip = 0.836(5), T c /m p = 0.224(1), and at k = 
0.19, P = 4.76, m^/m p = 0.708(7), T c /m p = 0.245(2). 
Using the lattice spacing estimated in the above, we find 
that at k = 0.190, f3 = 4.76, m,, = 578(2) MeV, at 
K = 0.160, P = 5.28, = 1991(7) MeV. Comparing 
the values of k, ft at simulation points in Ref. [37j with 
ours, we can roughly estimate the pseudo-scalar meson 
mass m v . Using the estimated lattice spacing, we can es- 
timate that Roberge- Weiss transition point temperature 
varies from 197 — 410 MeV in our simualtions. 

We consider the peak behaviour, reweighted distribu- 
tion and Binder cumulant of order parameter around 
the critical end point (3rw, At k = 0.190, 0.198, the 
three observables' behaviour show that transition at the 
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FIG. 7: Scaling behavior of susceptibilities of the imaginary part of the Polyakov loop according to the first order critical 
indexes (left panels) and to the 3D Ising critical indexes (right panels) . 
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FIG. 8: Reweighted distributions of the imaginary part of the Polyakov loop at k = 0.165 at the corresponding end points 
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end point is of first order which means the end point is a triple point. At k = 0.155,0.160, the peak behav- 
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FIG. 9: Reweighted distributions of the imaginary part of the Polyakov loop Im(L) at the corresponding end point /3rw, and 
P > Prw and /3 < Prw on each lattice spatial volume at k = 0.168. 



ior at the end point are more consistent with that of 
transition of a triple point than that of 3D Ising transi- 
tion behaviour. Similar observations can be observed at 
n = 0.170, 0.175, 0.180. 

At k — 0.165, 0.168, it becomes difficult to discern 
the peak behavior between 3D Ising transition class and 
triple point, however, when we look at the peak behavior 
at large lattice size, it is a reasonable conclusion that 
the behavior of RW transition line end points are of first 
order. This conclusion is enhanced by the reweighted 
distribution of order parameter. 

We also fit Eq. (TT21 to the calculated Binder cumulant 
data to extract the value of critical index v. At K = 
0.165, 0.168, v — 0.3661, 0.3594, respectively, and these 
values conform to first order transition. 

In Ref. [l3j], the locations of triple points are deter- 
mined. In Ref. [U, Q and Ref. Q , the simulations with 
staggered fermions show that phase diagram of two fla- 
vor and three flavor QCD at imaginary chemical potential 
fi = iirT are characterized by two tricritical points, re- 



spectively. Our simulations have no evidence that shows 
the existence of tricritical points separating second order 
region from the first order region. Considering these re- 
sults, our investigation requires further extensive numer- 
ical simulations which extend to a larger range of quark 
mass region. This work is under progress. 
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